#include "dyn_interpolation.h"

Interpolation::Interpolation(colvec _xValues, colvec _yValues, int _interPoints)
{
    // We set the interpolation variables
    setXValues(_xValues);
    setYValues(_yValues);
    setInterpolationPoints(_interPoints);

    // We set the back-end stuff variables.
    size = xValues.n_rows;
    correlation = 0;
    jsav = 0;
    dj = std::min(1,(int)pow((double)size,0.25));
}

double Interpolation::interpolate(double x)
{
    int index = correlation ? hunt(x) : locate(x);
    return rawInterpolation(index,x);
}

int Interpolation::locate(const double x)
{
    // We prepare the placeholders parameters.
    int lowerBound,upperBound,midPoint;

    // We check to see if we properly configured
    // the interpolation object.
    if (size < 2 || interPoints < 2 || interPoints > size)
        throw("Bad configuration of interpolation object.");

    // We determine if the sequence is monotically increasing
    // or monotonically decreasing.
    bool ascending = (xValues(size-1)>=xValues(0)); // True is increasing, false otherwise.

    // We initialize the bounds.
    lowerBound=0;
    upperBound=size-1;

    // We now look for the index having
    // the minimum \x-x_0\ value.
    while (upperBound - lowerBound > 1)
    {
        // Compute the middle point
        midPoint = (lowerBound+upperBound) >> 1;    // Bitwise operation ensures integer result.

        // If x_0 is bigger than x, set the midpoint to the lower
        // bound (for increasing sequences)
        if ((x >= xValues(midPoint)) == ascending)
            lowerBound = midPoint;
        else
            upperBound = midPoint;
    }

    correlation = abs(lowerBound-jsav) > dj ? 0 :1;
    jsav = lowerBound;

    return std::max(0,std::min(size-interPoints,lowerBound-((interPoints-2)<<1)));
}

int Interpolation::hunt(const double x)
{
    // We initalize the placeholder parameters.
    int lowerBound = jsav, upperBound, midPoint, stepsize=1;

    // We check to see if we properly setup the interpolation
    if(size < 2 || interPoints < 2 || interPoints > size)
        throw("Bad configuration of interpolation object.");

    // Increasing or decreasing sequence
    bool ascending = (xValues(size-1)>=xValues(0));

    if (lowerBound < 0 || lowerBound > size -1)
    {
        lowerBound = 0;
        upperBound = size-1;
    }
    // Hunt phase
    else
    {
        if ((x >= xValues(lowerBound)) == ascending) // Hunt up.
        {
            for (;;)
            {
                upperBound = lowerBound+stepsize;
                if (upperBound >= size-1){upperBound = size -1;} // We overshot the file.
                else if ((x < xValues(upperBound)) == ascending)
                    break;                                       // Found bracket.
                else
                {
                    // We are too low in the table. Let's look up and double
                    // the width of the area to search.
                    lowerBound = upperBound;
                    stepsize += stepsize;
                }
            }
        }
        else                                        // Hunt down.
        {
            upperBound = lowerBound;

            for (;;)
            {
                // Decrease the lower bound.
                lowerBound -= stepsize;
                if (lowerBound <= 0){lowerBound = 0; break;}    // Off end of table.
                else if ((x >= xValues(lowerBound)) == ascending)
                    break;                                      // Found bracket.
                else
                {
                    // We are too high in the table.
                    upperBound = lowerBound;
                    stepsize += stepsize;
                }
            }
        }
    }

    // Bisection phase.
    while (upperBound - lowerBound > 1)
    {
        // Compute the middle point.
        midPoint = (lowerBound+upperBound) >> 1;

        // Bisect.
        if ((x>=xValues(midPoint)) == ascending)
            lowerBound = midPoint;
        else
            upperBound = midPoint;
    }

    correlation = abs(lowerBound-jsav) > dj ? 0 : 1;
    jsav = lowerBound;
    return std::max(0,std::min(size-interPoints,lowerBound-((interPoints-2)>>1)));
}
